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ABSTRACT 

We perform numerical simulations of the growth of a Population III stellar system under photodissociating feedback. 
We start from cosmological initial conditions at z = 100, self-consistently following the formation of a minihalo at 
z — 15 and the subsequent collapse of its central gas to high densities. The simulations resolve scales as small as 
~1 AU, corresponding to gas densities of 10 16 cm -3 . Using sink particles to represent the growing protostars, we 
evolve the stellar system for the next 5000 yr. We find that this emerging stellar group accretes at an unusually low 
rate compared with minihalos which form at earlier times ( z = 20-30), or with lower baryonic angular momentum. 
The stars in this unusual system will likely reach masses ranging from < 1 Mq to ~5 Mq by the end of their 
main-sequence lifetimes, placing them in the mass range for which stars will undergo an asymptotic giant branch 
(AGB) phase. Based upon the simulation, we predict the rare existence of Population III stars that have survived to 
the present day and have been enriched by mass overflow from a previous AGB companion. 
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1. INTRODUCTION 

Following the emission of the cosmic microwave background, 
the universe entered a period referred to as the “dark ages,” 
when no luminous objects had yet formed. During this time, 
self-gravitating dark matter (DM) halos gradually grew in mass 
through the process of hierarchical merging. The dark ages 
ended when first stars, also known as Population III (Pop III), 
formed at z > 20 within DM minihalos of mass ~10 6 M 0 (e.g., 
Haiman et al. 1996; Tegmark et al. 1997; Yoshida et al. 2003). 

The typical mass of Pop III stars remains an open ques- 
tion that is crucial to our understanding of the evolution of the 
early universe (Bromm 2013). Their mass determines the rate 
at which they emitted ionizing radiation, and thus the extent to 
which they contributed to the reionization of the universe (e.g., 
Kitayama et al. 2004; Sokasian et al. 2004; Whalen et al. 2004; 
Alvarez et al. 2006; Johnson et al. 2007). In addition, the Pop III 
mass determines how much they contributed to the metallicity 
of the intergalactic medium (IGM; Madau et al. 2001; Mori 
et al. 2002; Bromm et al. 2003; Wada & Venkatesan 2003; 
Norman et al. 2004; Tornatore et al. 2007; Greif et al. 2007, 
2010; Wise & Abel 2008; Wise et al. 2012; Maio et al. 2011; 
recently reviewed in Karlsson et al. 2013). For instance. Pop III 
stars with main sequence masses in the range 40 Mq < 
M* < 140 M q or M* > 260 M 0 are expected to collapse di- 
rectly into black holes, therefore releasing virtually no metals 
into the IGM. On the other hand, stars with mass 140 M 0 < 
M* < 260 M q are predicted to explode as pair-instability su- 
pernovae (PISNe; Heger & Woosley 2002), thereby releasing 
the entirety of their metal content into the IGM and surrounding 
halos. We furthermore note recent work which has found that 
stellar rotation may significantly lower the PISN mass range 
(Chatzopoulos & Wheeler 2012; Yoon et al. 2012). Primordial 
stars within the range 15 M 0 < M* < 40 M 0 will end their 
lives as core-collapse SNe, or possibly hypernovae in the case 
of rapid rotation (e.g., Nomoto et al. 2003). Constraining the 
initial mass of Pop III stars is therefore central to understanding 


how the radiation and metallicity they released affected the for- 
mation of later stellar generations. 

Earlier work predicted that Pop III stars would form as 
single stars and grow to be very massive (>100 M 0 ; e.g., 
Abel et al. 2002; Bromm et al. 2002; Bromm & Loeb 2004; 
Yoshida et al. 2008). An analytical study by McKee & Tan 
(2008) found that even when accounting for radiative feedback, 
a protostar can grow to greater than 100 M 0 . Though a portion 
of the inflow toward the protostar will be reversed due to 
a growing ionization front (I-front), this I-front will expand 
preferentially in directions perpendicular to the protostellar disk, 
while accretion through the disk can continue unimpeded until 
much later times. 

In contrast to the above picture of single and massive Pop III 
stars, more recent work has shown that primordial gas will 
undergo fragmentation and develop into a disk within which 
a stellar multiple system will form (e.g., Clark et al. 2008, 
2011a). Such fragmentation is seen in simulations even when 
initialized on cosmological scales (e.g., Turk et al. 2009; Stacy 
et al. 2010). Furthermore, Pop III multiplicity occurs down to 
very small scales (~10 AU) and in the majority of minihalos 
(Clark et al. 2011b; Greif et al. 2011), even when accounting for 
the effects of feedback from protostellar accretion luminosity 
(Smith et al. 2011; Stacy et al. 2012). Though a number of the 
above simulations exhibited disk fragmentation, the resolution 
study by Turk et al. (2012) found that increasing the number 
of resolution elements per Jeans mass leads to variation in gas 
morphology and suppression of disk formation. They did not 
follow subsequent protostellar accretion, however, so whether 
fragmentation instabilities might later develop could not be 
determined. Latif et al. (2013) performed a similar study but 
followed the gas evolution beyond the formation of the first 
peak to find that self-gravitating disks with very rapid accretion 
rates ( 1 0~ 2 M 0 yr -1 ) will indeed develop. 

Several studies thus provide evidence that disk instability will 
develop after the first protostar arises in a primordial minihalo. 
In addition, a fraction of minihalos may also be subject to 
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earlier fragmentation in a pre-collapse phase. Greif et al. (2013) 
find that this will lead to secondary clumps in at least two 
out of nine minihalos even before a protostar forms. Though 
these recent studies generally imply a broader IMF for Pop III 
stars, the predicted IMFs remain top-heavy. The primordial 
stellar systems still exhibit rapid accretion rates compared to 
modern-day star formation (e.g., Stacy & Bromm 2013), and the 
most massive star in each system is still expected to eventually 
reach very high masses (>10 Mq). 

In this paper we present the first three-dimensional numerical 
simulation to follow the growth of a Pop III stellar system 
from cosmological initial conditions while also resolving nearly 
protostellar scales (~100 Rq) and accounting for the effects of 
photodissociating radiation. Recent work has found that 100 Rq 
(~1 AU) is approximately the maximum radius reached by a 
Pop III protostar during its pre-main sequence evolution (see 
e.g., Hosokawa et al. 2010; Smith et al. 2012). Resolving 
these small scales corresponds to evolving the gas up to a 
maximum density of 10 16 cnT 3 , at which point we continue 
the gas evolution for a further 5000 yr by employing sink 
particles to represent the growing protostars. At this high density 
the gas is quickly approaching complete optical thickness to 
continuum radiation (10 18 cm~ 3 ) and will not undergo further 
fragmentation on sub-sink scales (Yoshida et al. 2008). Our 
calculation therefore allows us to determine the true number 
of protostars that form within the central region of the host 
minihalo without missing any fragmentation due to lack of 
spatial resolution. 

While our numerical feedback model is also able to follow 
the effects of ionizing radiation from a growing protostar 
(e.g., Greif et al. 2009), we find that the particular protostellar 
system we simulate does not contain stars sufficiently massive 
to produce an Hn region. This is an unusual system with 
significant variation from the more typical rapidly accreting 
Pop III protostars studied in abovementioned work. Instead, the 
most massive star of the system considered here will most likely 
undergo an asymptotic giant branch (AGB) phase en route to 
becoming a white dwarf. 

The AGB phase, which occurs for stars with initial masses 
between ~0.8 and 8 M Q , plays an influential role in Galactic 
chemical evolution, so there is much interest in understanding 
the metal yields of low-metallicity AGB stars (e.g., Karakas 
2010; Campbell & Lattanzio 2008; Karakas & Lugaro 2010). 
Similarly, if a significant population of AGB stars existed at 
high redshift, this would have consequences for metal and 
dust production in the early universe. Observations of high- 
redshift (z > 6) galaxies and quasars indicate that significant 
amounts of dust had already formed at these early times, and the 
origin of such rapid dust production remains a subject of study 
(Bertoldi et al. 2003; Valiante et al. 2009; Cherchneff & Dwek 
2010; Gall et al. 2011). Though SNe were likely the major 
source of dust at these early times, Pop III AGB stars could 
have provided a significant contribution as well. This applies 
in particular to stars such as the three most massive ones we 
present from our simulation, predicted to reach 3-5 M 0 . These 
stars have sufficiently short main-sequence lifetimes (~10 8 yr) 
to undergo an AGB phase by z = 6. In contrast, the smaller 
1 Mq stars from our simulation are too long lived to provide 
any dust or metallicity contribution. The more massive AGB 
stars could also have significantly contributed to carbon and 
nitrogen production in the early universe, as well as s-process 
elements (Busso et al. 2001; Siess et al. 2002; Siess & Goriely 
2003). 


Though much study of Pop III stars to date has emphasized 
the high-mass end of the Pop III IMF, we note that even 
some early studies predicted that typical Pop III stellar masses 
might be quite low, <1 Mq. For instance, Kashlinsky & Rees 
(1983) emphasized the importance of angular momentum in 
determining the mass of Pop III stars, predicting that rotational 
effects would cause the primordial gas clouds to collapse into a 
dense disk. Only after the disk cooled to ~1000 K through 
FF line emission could fragmentation occur. Nakamura & 
Umemura (2001) predicted that fragmentation of primordial 
filaments would lead to a bimodal IMF, with a peak at ~ 1 Mq 
as well as at 100 Mq. 

Here, we further explore the possible parameter space for 
Pop III star formation by modeling with high accuracy the 
growth of an unusually low-mass primordial system. Such 
cases are expected to be rare, since most Pop III systems have 
been found to contain one or more high-mass stars that likely 
dominated the overall Pop III impact on the IGM and later 
stellar generations. However, such low-mass stars as found in 
our simulation were potential survivors to the present-day, and 
may in principle be discovered within the Milky Way halo or 
nearby dwarf galaxies. In Section 2 we describe our numerical 
methodology, in Section 3 we discuss our protostellar evolution 
model, and in Section 4 we present our results. We discuss 
the impact of a global Lyman-Werner (LW) background in 
Section 5, and we conclude in Section 6. 

2. NUMERICAL METHODOLOGY 
2.1. Initial Setup 

We carry out our investigation using gadget-2, a widely 
tested three-dimensional /V-body and smoothed particle hydro- 
dynamics (SPH) code (Springel 2005). We begin with a 200 kpc 
(comoving) box containing 128 3 SPH gas particles and the same 
number ofDM particles. The simulation is initialized at z = 100. 
Positions and velocities are assigned to the particles in accor- 
dance with a ACDM cosmology with £2 A = 0.7, = 0.3, 

£2 b = 0.04, os = 0.9, and h = 0.7. The gas and DM evolu- 
tion is followed until the first minihalo forms and its central gas 
density reaches 10 4 cm -3 . 

Once the site of the first minihalo is determined, the simula- 
tion is performed at higher resolution, starting again at z — 100. 
The increased resolution is attained using a hierarchical zoom- 
in procedure (e.g., Navarro & White 1994; Tormen et al. 1997; 
Gao et al. 2005) in which four nested refinement boxes are 
placed within the cosmological box, centered on the site where 
the minihalo will eventually form. Within each refinement level, 
each particle from the lower level is replaced with eight “child” 
particles, so that in the most refined region the parent particle is 
replaced by 4096 child particles. The four refinement levels have 
lengths of 40, 35, 30, and 20 h~ x kpc (comoving), so that the 
most highly refined level encompasses all the mass that will later 
be incorporated into the minihalo. The most refined SPH parti- 
cles have mass msph = 5 x 1 0 3 M 0 , and the resolution mass 
of the refined simulation is M res ~ 1 .5/V ne j g h/nspH < 0.3 Mq, 
where A ne i g h — 40 is the typical number of particles in the SPH 
smoothing kernel (e.g., Bate & Burkert 1997). 

2.2. Cut-out Technique and Particle Splitting 

To increase computational efficiency, once the gas has reached 
densities of 10 12 cm -3 we employ a “cut-out” technique in 
which all gas and DM beyond 3 pc from the densest gas 
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particle is removed. At this point the central star-forming gas 
is gravitationally bound and under minimal influence from the 
mass of the outer minihalo and more distant regions of the 
cosmological box. The total mass of the cut-out is 3500 M 0 , 
and the minimum density is ~10 2 cm~ 3 (see, e.g., Stacy et al. 
2012 for further details). 

At the same time that we cut out the central 3 pc of 
the cosmological box, we also further increase the particle 
resolution so that collapse to densities of 10 16 cm -3 can be 
properly followed. We thus replace each SPH particle with 
eight child particles, each of which is placed randomly within 
the smoothing kernel of the parent particle. The mass of the 
parent particle is then evenly divided amongst the child particles. 
Each of these particles inherits the same chemical abundances, 
velocity, and entropy as the parent particle (see, e.g., Bromm 
& Loeb 2003; Clark et al. 2011b). This ensures conservation 
of mass, internal energy, and linear momentum. After this 
process, each SPH particle in the new cut-out simulation has 
a mass of m sp h = 6 x 1 0 4 M 0 , and the new resolution mass 
is M ies — 0.03 M 0 . This final M res is close to the mass of the 
pressure-supported atomic core which develops once the opacity 
limit is reached (Yoshida et al. 2008), defining the point at which 
the protostar has first formed. 

2.3. Chemistry, Heating, and Cooling 

We utilize the same chemistry and thermal network as de- 
scribed in detail by Greif et al. (2009) and used in Stacy et al. 
(2012). In short, the code follows the abundance evolution of 
H, H + , H _ , H2, Ht, He, He + , He ++ , and e~ , as well as the three 
deuterium species D, D + , and HD. All relevant cooling mecha- 
nisms, including H2 collisions with H and He as well as other 
H2 molecules, are included. The thermal network also accounts 
for cooling through H2 collisions with protons and electrons, 
H and He collisional excitation and ionization, recombination, 
bremsstrahlung, and inverse Compton scattering. 

Further H2 processes must also be included to properly model 
gas evolution to high densities. For instance, the chemistry and 
thermal network includes three-body H2 formation and the con- 
comitant H2 formation heating, which become important at 
n > 10 8 cm -3 . Furthermore, when n > 10 9 cm -3 , cooling 
through H2 ro-vibrational lines becomes less effective as these 
lines grow optically thick. This is accounted for using an escape 
probability formalism together with the Sobolev approximation 
(see Yoshida et al. 2006; Greif et al. 2011 for further details). 
Simple fitting formulae are also available for estimating opti- 
cally thick H2 rates (e.g., Ripamonti & Abel 2004). However, 
Hirano & Yoshida (2013) find that fitting formulae can over- 
estimate the cooling rate, such as in cases when the gas has 
a substantial degree of rotation. This can lead to significant 
differences in gas evolution, such as accelerated gravitational 
collapse, when using fitting formulae as opposed to the more 
accurate Sobolev method. 

The most important new process utilized in the thermal 
network is H2 collision-induced emission (CIE) cooling, which 
becomes significant when n > 10 14 cm -3 (Frommhold 1994). 
As described in Greif et al. (2011), the reduction of the CIE 
cooling rate due to the effects of continuum opacity is handled 
through the following prescription (Ripamonti et al. 2002; 
Ripamonti & Abel 2004): 

( 1 - e~ Tcm \ 

Acie, thick = AciE.thin min I , 1 I > (1) 

V t cie ) 


where 

r cm = ( a > ( 2) 

\7 x 10 13 cm 7 

Acie, thin is the CIE cooling rate in the optically thin limit, and 
Acie, thick that in optically thick conditions. Hirano & Yoshida 
(2013) compared gas evolution when continuum opacity effects 
are calculated using the above fitting formula, and when they 
are instead estimated using three-dimensional ray tracing. They 
find that when the fitting formula is used the gas collapses to 
~10 17 cm -3 only slightly faster, by ~ 1 yr. The differences in the 
thermal evolution are also minimal between the two methods. 
We thus expect the fitting formula above to be sufficiently 
accurate. However, while Hirano & Yoshida (2013) modeled 
runaway gas collapse, it is possible that the differences would 
be more substantial when considering longer-term evolution of 
a disk, and this will be further examined in future work. 

2.4. Sink Particle Method 

When an SPH particle reaches a density of n max = 10 16 cm -3 , 
it is converted to a sink particle (e.g.. Bate et al. 1995; Bromm 
et al. 2002; Martel et al. 2006; see also Stacy et al. 2012 for 
further details on the employed sink particle method.) The sink 
then grows in mass by accreting surrounding particles within 
its accretion radius, which we set equal to the resolution length 
such that r acc = L res ~ 1.0 AU. 

The sink accretes a gas particle within r acc as long as the 
particle is not rotationally supported against infall onto the 
sink. This is determined by checking that the particle satisfies 
./SPH < 7 cent? where ;' SP h = v mt d is the specific angular 
momentum of the gas particle, y cent = y/G M sm ^r. dcc the level 
required for centrifugal support, and u rot and d are the rotational 
velocity and distance of the particle relative to the sink. Particles 
that satisfy these criteria are removed from the simulation, and 
their mass is added to that of the sink. When the sink first 
forms it immediately accretes most of the particles within its 
smoothing length, so its initial mass is near the resolution mass 
of the simulation, M res ~ 3x 1 (T 2 M 0 . Its position and velocity 
are set to the average of that of the accreted particles. 

These same accretion criteria are additionally used to deter- 
mine whether two sinks may be merged. However, we note that 
modifications to the sink merging algorithm can significantly al- 
ter the sink accretion history (Greif et al. 2011). Recent work by 
Greif et al. (2012) has resolved sub-protostellar scales (0.05 R Q ) 
of primordial star-forming gas, tracking the merger rate of pro- 
tostars by following their interactions without using sinks. They 
find that approximately half of the secondary protostars formed 
will indeed migrate toward and merge with the initial protostar. 
Our merging algorithm leads to a similar fraction of secondary 
sinks merging with the main sink. After the first sink arises, six 
secondaries later form, but two of them eventually merge with 
the initial sink. This roughly agrees with what Greif et al. (2012) 
find to occur on sub-sink scales. 

2.5. Ray-tracing Scheme 

Once the first sink particle forms, it represents a newly formed 
protostar and is used as the point source for modeling the effects 
of LW radiation emanating from the protostar. We use the same 
scheme as described in Stacy et al. (2012). Briefly, our ray- 
tracing module generates a spherical grid, consisting of ~10 5 
rays and 200 radial bins, centered around the first sink. The 
minimum radius is set equal to the distance between the sink 
and the nearest neighboring SPH particle, and the grid is updated 
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each time ray tracing is performed. Most particles within r acc 
from the sink are accreted, so the minimum radius is usually 
close to 1.0 AU. The bins are logarithmically spaced from 
the minimum distance to 3 pc, the size of the cut-out region. 
Each SPH particle within a bin then contributes its density and 
chemical abundances, proportional to its density squared, to the 
average values assigned to the bin. 

Also as in Stacy et al. (2012), we then use the ray-tracing 
scheme to determine the H 2 column density /Vh 2 , and from 
this we determine the shielding factor /shield with the fitting for- 
mula from Draine & Bertoldi (1996). We note a recent update 
to the /shield prescription by Wolcott-Green et al. (2011) and 
Wolcott-Green & Haiman (201 1), but do not expect this to sig- 
nificantly affect the results for our particular case. We combine 
this with a protostellar evolution model (see Section 3) in which 
we assume a blackbody spectrum with an effective temperature 
Teff, as specified by that model. We then determine the approx- 
imate LW radiation flux 7/ w, in units of erg s 1 cm~ 2 Hz -1 , 
at hv — 12.87 eV (Abel et al. 1997). This finally allows for a 
determination of the H 2 dissociation rate, 

k H2 = 1.1 X 10 8 /shield F lw s’ 1 , (3) 

to be included in our chemical network. 

3. PROTOSTELLAR EVOLUTION MODEL 
3.1. Luminosity and Temperature Evolution 

Our ray-tracing algorithm requires an input of protostellar 
effective temperature T e ff and luminosity L.,. We calculate L* as 
the sum of L acc , the accretion luminosity, and L; nt , the luminosity 
originating from the stellar interior and finally emitted from the 
photosphere of the protostar. We write L* as 

GM*M 

L* = L acc + Ei n t = ot — I- L m t, (4) 

where M* is the protostellar mass, M the accretion rate, and 
R* the protostellar radius (cf. Prialnik & Livio 1985; Hartmann 
et al. 1997). We take M* to be the mass of the sink, and M to 
be the accretion rate onto the sink, measured by averaging the 
total mass growth of the sink over the previous 10 yr. If the 
measured sink accretion rate yields M ~ 0, we simply assume 
the protostar is described by L* = L; nt . 

We define a as the fraction of thermal energy from accretion 
that is added to the stellar interior. Cold disk accretion would 
thus be described by a = 0, while for hot spherical accretion 
a — 1 . For the main sink particle, we determine a by measuring 
the /sph of each particle accreted by the sink within the last 
10 yr, as well as each particle currently within 10 AU from 
the sink. This allows us to find the percentage of nearby and 
recently accreted particles that have low angular momentum 
Osph < 0.5 jcent) versus high angular momentum O'sph > 
0.5 ,/ C ent)- As described in Section 2.5, /sph = t> rot d is the angular 
momentum of the gas particle and / cen t = \/G is 

the angular momentum required for centrifugal support against 
infall onto the sink. We take a to be the fraction of particles 
with /sph < 0.5 / C ent, such that a = 1 if the near-sink gas is 
dominated by radial instead of rotational motion. 

L; nt will vary with the mass of the protostar. When the 
protostar first forms and has low mass, we assume it is 
on the Hayashi track of the Hertzsprung-Russell diagram. 
We approximate this by holding the protostar at an effective 


temperature of Tkay = 4500 K while its luminosity may vary. 
This yields a “Hayashi track” luminosity of 

LHay = 4lT R-OssT^y. (5) 

The precise value of 7nay will vary from ~3000 to 5000 K de- 
pending upon stellar mass and opacity. Because the protostar 
we consider is metal-free, the opacity of the protostellar atmo- 
sphere will differ from that of a Pop I or Pop II star, and the 
resulting 7 Hay will tend to be marginally higher for Pop III stars 
(e.g., Stahler et al. 1986a). We thus choose ?Hay to be in the 
upper end of this range and set it to 4500 K. The uncertain value 
of initial T e g- (from 3000 to 5000 K) corresponds to a Lh ;iv that 
may vary by a factor of eight. However, a 7/n even in the upper 
end of this range will not lead to significant feedback until after 
the protostar leaves the Hayashi track. Our simulation results 
are thus not sensitive to the choice of initial L e ff. 

If the protostar grows sufficiently massive, we assume it 
eventually transitions to the Henyey track (Henyey et al. 1955, 
see also Hansen et al. 2004), and will gradually contract down 
to the main-sequence radius and commence hydrogen burning. 
However, the protostellar system in our simulation exhibits 
unusually low mass, and thus has much longer evolutionary 
timescales than the more common high-mass Pop III stellar 
systems. We thus do not follow sufficiently long timescales 
for the protostars to reach these later stages, so we do not 
discuss them here. Our protostellar model at these early times 
predicts a typical L* and 7/n of ~100 L Q and 4500 K. This 
luminosity can roughly double during periods of rapid accretion 
due to the contribution from L ac c , while the corresponding 7/n 
will increase by a few hundred kelvins. At such low L^ff, the 
fraction of luminosity in the LW band is negligible. Combined 
with significant H 2 shielding in the disk, LW dissociation is 
unimportant at these early times. 

3.2. Radial Evolution 

As discussed in previous one-dimensional studies (Stahler 
et al. 1986b; Omukai & Palla 2003; Hosokawa et al. 2010), a 
growing protostar initially undergoes an “adiabatic accretion” 
phase, where R* grows with mass. This continues approximately 
while f acc < fKH, where 


Lch 


CM; 

R*L* 


is the Kelvin-Helmholtz (KH) timescale and 


t 


acc 


M* 

M 


( 6 ) 


(7) 


is the accretion timescale. The protostar will later begin KH 
contraction approximately when f acc > r K H- 

Hosokawa et al. (20 1 0) modeled the evolution of a primordial 
protostars growing at M = 10 -3 M 0 yr -1 . They found that, 
particularly during the “adiabatic” phase, accretion through 
a geometrically thin disk will lead to smaller protostellar 
radii than spherically symmetric accretion. The true accretion 
geometry is likely somewhere in between the idealized “disk” 
and “spherical” cases, with the infall beginning as nearly 
spherical and growing gradually more disk-like over time. Given 
the spherical accretion case described in Omukai & Palla (2003; 
see also Stahler et al. 1986b), the radial evolution during the 
adiabatic accretion phase can be described by the following 
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/ M* \ 1/: 3 / M \ 1/3 

' 49 *U) (sj ' (8) 

where M gd — 4.4 x 10” 3 M 0 yr” 1 , the fiducial accretion rate 
used in the abovementioned studies (see also Stacy et al. 2010). 

If a protostar transitions from spherical to disk accretion, the 
radial evolution during the “adiabatic accretion” phase will be 
significantly different from the purely spherical case. For the 
range of accretion rates studied in Hosokawa et al. (2010), after 
transitioning to disk accretion, the radius rapidly declines due to 
the decrease in entropy brought to the stellar interior. For pure 
disk accretion, this decline can be described by 

/M*V°' 63 

^ sk ^°Uo) ’ W 

where Ro and Mo are the protostellar mass and radius at the point 
of transition from spherical to disk accretion. In our case, the 
protostar’s accretion is disky even from initial sink formation, 
so we simply use the initial sink mass to set Mq = 0.045 M Q . 
Following the disk accretion model of Hosokawa et al. (2010; 
their Figure 4), we approximate Rq — 1 ,7M* _ l/3) = 5 Rq. 

We subsequently set our radius in between the spherical and 
disk cases such that 

R / — sphere (1 C*0 Rl, disk- (10) 

If the radial decline is unphysically rapid ( R < — R/t KH ), we 
limit the rate of contraction to be R = — R/t-^a- The radial 
decline will continue until deuterium burning begins in the 
stellar interior, thus increasing the average entropy within the 
star, which occurs when 7i nt > 2 x 10 6 K. After this point 
the protostar again undergoes a roughly adiabatic expansion. 
When a maximum radius is reached (see Hosokawa et al. 2010), 
KH contraction to the main sequence will begin. However, our 
simulations do not follow the protostellar growth to these later 
stages because we are examining an atypical low-mass case 
in which the protostellar evolutionary timescales are very long 
compared to those of more massive protostars. 

In our model R * expands to nearly 30 R Q over the first few 
hundred years during a period of more spherical-type accretion. 
After a transition to more disk-type accretion, R t gradually 
declines as R — —R/ r KH down to ~ 1 5 R Q by the end of the 
simulation. 

4. RESULTS 

4.1. Initial Minihalo Collapse 

The chemical and thermal evolution of the central minihalo 
gas up to the formation of the first sink particle is depicted by 
the red lines in Figure 1. The minihalo is in place by z = 15, 
and the subsequent evolution is similar to that of previous 
work (e.g., Yoshida et al. 2006; Greif et al. 2011). The gas is 
heated through adiabatic compression as it approaches densities 
of 10 s cnT 3 . After this point, three-body reactions rapidly 
increase the H 2 fraction, such that the correspondingly enhanced 
cooling rate is similar to the combination of the H 2 formation 
heating rate and adiabatic heating rate due to collapse, yielding a 
roughly isothermal evolution. The gas becomes fully molecular 
by densities of 10 10 cm" ', forming a ~1 M 0 molecular core. 
Above these densities the H 2 cooling is no longer optically thin 
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and the gas gradually heats again to ~2000 K by n = 10 I6 cm” 3 , 
at which point the first sink particle forms. 

The velocity structure of the central 10,000 AU is shown 
in Figure 2. Within each logarithmically spaced radial bin, we 
take u ra d and v ro t as the mass-weighted average of the individual 
particle velocities within each bin. Both the radial and rotational 
gas velocity u ra d and i> rot are on the order of half of vg, where 
riff is the free-fall velocity based upon the enclosed mass at the 
given radius. The gas has a substantial amount of rotational 
support such that ?; rl)t dominates over u ra d and is approximately 
half of the Keplerian velocity i)K ep - 

In a similar fashion, we measure the turbulent Mach number, 
Mfurb s over the same range of radial bins according to 

M mb c l = _ "rad) 2 , (ID 

i 

where c s is the sound speed of the radial bin, m, is the mass of a 
gas particle contributing to the bin, and M is the total gas mass 
of the bin. The central 10,000 AU of gas are characterized by 
subsonic and nearly sonic turbulence. 

4.2. Comparison with Other Minihalos 
4.2.1. Global Minihalo Characteristics 

As will be discussed in Section 4.3, our simulated Pop III 
system has an unusually low accretion rate. We here examine 
whether this is due to the characteristics of its host minihalo. 
We first measure the evolution of the virial mass of the 
minihalo considered here, and compare to four other minihalos 
taken from the cosmological simulations presented in Greif 
et al. (2012) and Stacy & Bromm (2013), where we have 
chosen the minihalos hosting the most rapidly and most slowly 
accreting stellar systems from each of those two studies. We 
determine which particles reside in the halo by first locating 
the simulation region’s densest gas particle, or hydrodynamic 
mesh element in the case of theAREPO simulations. Making 
the simple assumption that this point marks the center of the 
halo, the extent of the halo was determined by finding the 
surrounding spherical region in which the average DM density 
is 200 pb, where pb — 2.5 x 10” 30 (1 + z) 3 g cm -3 is the redshift- 
dependent background density. 

The minihalo of our simulation has the minimum necessary 
mass of Mhaio 10 6 Mq before gas condensation and H 2 
cooling begins, and this is very close to the mass of other Pop III 
star-forming halos (Figure 3). However, it does not reach this 
minimum M^io until the relatively late time of z = 15. By 
this redshift, Hubble expansion has reduced the density of the 
background universe, leading to a slower accretion rate from the 
surrounding cosmic web. In particular, the average DM accretion 
rate Mdm of our minihalo is 3 x 1 0” 3 M 0 yr” 1 , as compared 
to the more typical rate of 10” 2 M 0 yr” 1 for minihalos which 
collapsed at z ~ 30. Employing this same redshift range and 
Mdm = 5 x 10 s Mq, we find that our measured Mdm values 
are in excellent agreement with analytical fits determined from 
other numerical simulations, e.g., 

Mdm ^ 35(1 + z) 2 ' 2 M,' 2 07 M 0 yr” 1 , (12) 

where M\i = Mdm/10 12 Mq (Genel et al. 2008). Though 
the above estimate was originally derived from simulations 
of much larger-mass halos, it also applies accurately to our 
minihalos. Genel et al. (2008) furthermore find that the above 
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Figure 1. Physical state of the minihalo gas 3000 yr after the formation of the first sink particle, (a) Temperature vs. number density n. (b) n vs. distance r from the 
main sink particle, (c) H 2 fraction /h2 vs. n. (d) H 2 fraction vs. distance r. Red line shows the radially averaged values of the same quantities just prior to the formation 
of the first sink, as measured from the densest gas particle. Note the warm phase of gas at n > 10 7 cm -3 , where the gas has been warmed through gravitational heating 
provided by the main sink. 

(A color version of this figure is available in the online journal.) 



Figure 2. Velocity profile of the gas just prior to initial sink formation. Solid 
line is the magnitude of the radial infall velocity u ra( j, averaged within a range 
of logarithmic radial bins. Dashed line is the average rotational velocity u ro t 
within these same bins. Dash-dot line is the free-fall velocity Uff based upon the 
enclosed mass at each radius. Dotted blue line is the turbulent Mach number 
Mturb, which corresponds to the right-hand y-axis scale. The values of u ra d are 
low compared to i>ff, and the magnitude of u ro t dominates over u ra d over all 
distances shown. 

(A color version of this figure is available in the online journal.) 


expression agrees well with the general predictions of extended 
Press-Schechter (EPS) theory (Press & Schechter 1974; Mo 
& White 1996; Lacey & Cole 1993, see also Neistein et al. 
2006). As numerically confirmed in Gao et al. (2005; e.g., their 
Figure 1), the Press-Schechter formalism indeed accurately 
predicts the accretion rate of minihalos with Mhaio ~ 10 6 M 0 . 
They found an accretion rate of ~2.5 x 10~ 2 M 0 yr -1 for 
their z = 50 minihalo when it was collecting the bulk of its 
mass, which in turn corresponds well with the above equation’s 
prediction of ~2 x 10~ 2 M 0 yr -1 . This further confirms the 
expected larger growth rates of these rare-peak and high-redshift 
minihalos (e.g., Reed et al. 2005). 

The gas mass within the minihalo, which we take simply as 
all gas within the halo virial radius with n > 1 cm -3 , grows at 
average rates which range from M gas = 1 0 4 M 0 yr -1 for our 
minihalo to 1 0 3 M 0 yr -1 for the earlier-collapsing minihalos. 
This may be compared to the growth rate found in previous 
studies, e.g., 

M m ~ 6.6 (1 + z) 125 M; 2 l 5 /o.i 65 M q yr” 1 , (13) 

where /o.i 65 is the baryonic fraction in the halos in units of 
the cosmological value /b = Ob/^m = 0.165 (e.g., Dekel 
et al. 2009; see also Faucher-Giguere et al. 2011). Assuming 
/o !65 = 1 we find that, over the redshift range z — 15-30, M gas 
ranges from 2 to 7 x 1 0 -4 M 0 yr -1 . Thus, for a given DM halo 
mass (e.g., Mh a i 0 = 10 6 M 0 ), the rate at which both DM and gas 
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Figure 3. Left: evolution of the DM mass of various simulated halos over time, as measured in terms of the age of the universe. Right: concurrent evolution of the 
baryonic mass within the halos. The red line represents the z = 15 minihalo of our simulation. The black dotted line represents the Stacy & Bromm (2013) minihalo 
which had the lowest overall stellar accretion rate, while the black dashed line represents that which had the greatest accretion rate. Green dotted and dashed lines 
additionally show the growth of the two minihalos from Greif et al. (2012) which had the lowest and highest accretion rates, respectively. 

(A color version of this figure is available in the online journal.) 


will fall into the gravitational potential well will vary by nearly 
an order of magnitude over this redshift range. This corresponds 
to the nearly factor of ten reduction in pb as the universe expands 
over redshifts z = 30 to z = 15. 

The variation in collapse redshifts seen in these simulations 
stems from differences in sizes of the cosmological box that was 
used. The Greif et al. (2012) and Stacy & Bromm (2013) simula- 
tions employed box sizes of 500 kpc and 1 .4 Mpc (comoving), 
and were thus able to capture larger “va peak” fluctuations, 
where a is the standard deviation in the Gaussian random field 
of primordial density fluctuations. ACDM theory predicts that 
at z ~ 30, a 10 6 M 0 halo corresponds approximately to a 3 a 
peak (e.g., Loeb 2010; Bromm 2013). Our smaller box size of 
200 kpc (comoving) captures only ler-2a peaks, which instead 
corresponds to a later collapse redshift of z ~ 15 for a 10 6 M 0 
halo. 

4.2.2. Star-forming Core Characteristics 

Along with slower overall halo gas accretion in the z = 
15 minihalo, we find a reduced gas accretion rate within the 
central parts of the minihalo as well, which may be roughly 
estimated through the gas soundspeed c s (Figure 4). Although 
we find good agreement in c s beyond 10 5 AU (~0.5 pc), there is 
divergence in c s in the inner regions. In our simulation the sound 
speed of the gas within 10 4 AU is ~3 km s' 1 . This is slightly 
less than that seen in some of the comparison halos, ~4 km s' 1 
(Figure 4). Making the simple estimate that the Jeans mass is 
infalling at the free-fall rate, we can scale the accretion rate with 
soundspeed as M ~ c 3 / G, similar to the Shu (1977) similarity 
solution for collapse of a singular isothermal sphere. From this 
we may predict accretion rates ranging from ~6 x 1 0 3 M 0 yr' 1 
for the z = 15 minihalo to ~1.5 x 10' 2 M 0 yr' 1 for the higher- 
redshift halos. 

Other properties of the gas, such as the density profile 
and FU fraction (lower panels of Figure 4), show interesting 
variation between minihalos, as well. At distances greater than 
100 AU, the H 2 fraction can vary by approximately an order 
of magnitude, though in each star-forming cloud the gas is 
fully molecular within the central 100 AU. We do not find 
an exact correlation between FU fraction and c s . However, our 


Z = 1 5 minihalo generally has the lowest Hi fraction at distances 
greater than 1000 AU as well as nearly the lowest temperatures 
at all distances. The gas density is unusually small as well. 
Though for each star-forming core the density roughly follows 
a p oc r' 2 profile, beyond 20 AU the z = 15 minihalo profile is 
normalized to smaller values than the others. At some radii the 
gas density is an order of magnitude lower than the minihalo 
with the highest density. 

As shown in Figures 4 and 5, we may also approximate the 
spherical accretion rate M sp h that results from the density and 
radial velocity profiles, which is appropriate when considering 
gas within the minihalo where the density profile is more 
spherically symmetric. We estimate A7 sph within a range of radial 
bins as 

M sp h = 4itr 2 pv rdd , (14) 

where r is the distance as measured from the densest gas particle, 
and i> ra d is the average radial velocity of gas within the radial bin. 
As is apparent in Figure 4, between 100 AU and 10 5 AU (0.5 pc) 
the gas within our minihalo typically infalls at v rdd ~ 1 km s' 1 , 
several times smaller than i> ra d ~ 4 km s' 1 as seen in the most 
rapidly growing halos. Together with the gas density this yields 
a gas accretion rate between 100 and 10 5 AU that ranges from 
~6 x 10' 4 M q yr' 1 to ~3 x 10' 3 M 0 yr' 1 , approximately an 
order of magnitude less than the accretion rates found in the 
high-redshift halos. For all minihalos, M sp h is also substantially 
lower than the accretion rate estimated from c s , indicating that 
angular momentum support slows gas infall. We furthermore 
note that the halo which formed at the second-most recent 
redshift (green dotted lines in Figures 3-5) has DM and gas 
accretion rates and M sv h values which are intermediate between 
the z = 15 minihalo and the highest-redshift halos. 

These results may indicate a correlation between collapse 
redshift and rate of mass infall even on small scales. However, 
Gao et al. (2007) found that while gas within higher-redshift ha- 
los will indeed reach protostellar densities in shorter timescales, 
there was no subsequent correlation between formation red- 
shift and instantaneous accretion rate at the densities they re- 
solved (10 10 cm' 3 ). O’Shea & Norman (2007) even find that 
Pop III star-forming regions which form later have higher max- 
imum accretion rates, the opposite trend to ours. They attribute 
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Figure 4. Radial profiles of various gas properties with respect to the highest-density particle, just prior to the initial protostar or sink particle formation. In all panels, 
dashed and dotted lines represent minihalos from simulations presented in Stacy & Bromm (2013) and Greif et al. (2012), as described in the previous figure. Upper 
left: number density vs. radius. Upper right: H 2 fraction vs. radius. Lower left: profile of sound speed c s . Lower right: profile of u ra d- The combination of smaller 
sound speed, as well as lower |u ra dl and number density, lead to unusually low accretion rates within our z = 15 minihalo as compared with those from other studies. 
(A color version of this figure is available in the online journal.) 


this to the increasing virial temperature of halos with redshift, 
7’vi, oc M v / r (1 + z). In their chain of reasoning, this would lead 
to warmer overall gas temperatures as the minihalo first col- 
lapses, yielding more rapid H 2 formation rates. The higher H 2 
fraction in turn would lead to cooler gas in the cores of the ha- 
los, and thereby slower accretion rates. In Figure 4, however, we 
find a counterexample to this. A relatively high-redshift mini- 
halo with an enhanced H 2 fraction (black dashed line) in fact 
has the warmest temperatures in both the core and outer region, 
and in some regions beyond 100 AU, it also has the highest 
accretion rates. 

Gao et al. (2007), on the other hand, emphasize the importance 
of angular momentum. They find that as gas collapses its 
properties become independent of the global properties of the 
halo, and that more disk-like and rotationally supported inner 
star-forming clouds will have lower accretion rates. In the 
bottom panel of Figure 5 we compare the levels of rotational 
support throughout the halos, defined as / rot = u rot /i>Kep, where 
fKep = (GM enc /r) 1/2 . The z = 15 minihalo has f m \ ~ 0.6 
over several orders of magnitude in distance. The other two 
more slowly accreting systems (green and black dotted lines 
in Figure 5) also maintain high/ rot between 10 and 10 6 AU, 
with/ rot consistently greater than 0.6. In contrast, the two most 
rapidly accreting systems (green and black dashed lines) have 
extended regions where f mt falls below 0.5. Thus, similar to the 
conclusions of Gao et al. (2007), angular momentum structure, 
combined with that of temperature and density, plays a key role 
in the rate of gas infall. 


Despite variation in the level of rotational support, the specific 
angular momentum profile of the central gas for each minihalo 
we examine, shown in Figure 6, all follow a similar j tot oc M enc 
power law. However, it is interesting to note that the profiles of 
the two most rapidly accreting halos (green and black dashed 
lines in Figure 6) are normalized at nearly a factor of two lower 
than the minihalo of our simulation. The greater rates at which 
the central gas flows inward (Figure 4), as well as the reduced 
rotational support (Figure 5, bottom panel), indicates that in 
these halos i> ra d dominates over i> rot , leading to higher accretion 
rates at the point of protostellar formation. 

As expected, the mass enclosed within a given radius is also 
larger for those minihalos with less rotational support. The 
resulting ratio of enclosed mass M enc to Bonnor-Ebert mass 
Mbe is larger as well (Figure 7). We estimate Mbe as 

/ T \ 3/2 / n \-i/2 

M BE - 1000 M 0 (— ; A . (15) 

° \ 200 K / V 10 4 cm~ 3 / 

Mbe is similar to the Jeans mass, and a ratio of M enc to Mbe that 
is close to unity roughly indicates gravitational instability of 
the gas. 

We note that de Souza et al. (2013) also emphasize the 
influence of rotation on the Pop III IMF, finding from their 
cosmological simulation that the spin distribution of gas within 
minihalos evolves with redshift. They employed the model from 
McKee & Tan (2008) which used semi-analytic methods to 
find a relation between gas rotational support, effectiveness of 
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Figure 5. Top: estimated spherical accretion rate M sp h over a range of distances 
from densest gas particle, just prior to protostar or sink particle formation. Line 
styles have the same meaning as in previous figures. Note that the combination 
of reduced density and |u ra dl within the inner few 10 s AU (~0.5 pc) leads to 
lower accretion rates than typically found in primordial gas. Bottom: level of 
rotational support / rot = c r ot/' ; Kcp throughout the gas at this same time. Note 
that the two most rapidly accreting systems (green and black dashed lines) are 
also the only two systems with extended regions of very low rotational support 
(/rot < 0.5). 

(A color version of this figure is available in the online journal.) 

protostellar feedback, and ultimate protostellar mass. Assuming 
one star per minihalo, de Souza et al. (2013) determine that the 
Pop III IMF should evolve to have lower peak masses at lower 
redshift. However, they point out that correlating spin with a 
Pop III IMF will be complicated by further feedback effects and 
Pop III multiplicity. 

Our results concerning the relation between angular momen- 
tum and protostellar accretion rate agree with the those of Gao 
et al. (2007). The implication that the IMF may shift to lower 
mass with lower redshift also shows a rough agreement with de 
Souza et al. (2013), but for different physical reasons. However, 
neither these authors nor O’Shea & Norman (2007) simulated 
the evolving accretion rate of the protostellar cloud after the 
first protostar formed. Our study allows for this through the sink 
particle method. As will be further discussed in Section 4.3, we 
find that the slower and more rotationally dominated gas infall 
within the z = 15 minihalo later leads to reduced accretion rates 
on to the evolving stellar system as well. In Figure 8 we compare 
the total sink mass after 5000 yr and the total angular momen- 



Figure 6. Angular momentum profiles from various simulated minihalos, taken 
at the point of first sink or protostar formation. The solid red line is taken 
from the simulation discussed here. The dashed black line is the most rapidly 
accreting cluster from Stacy & Bromm (2013). The dotted black line is the most 
slowly accreting cluster from Stacy & Bromm (2013). The green dotted and 
dashed lines are taken from Greif et al. (2012) minihalos which had the lowest 
and highest accretion rates, respectively. 

(A color version of this figure is available in the online journal.) 

turn within the central 200 M G just prior to sink formation. We 
use the minihalo of this work as well as the suite of ten miniha- 
los presented in Stacy & Bromm (2013). Some anti-correlation 
is apparent, particularly in that the z = 15 minihalo has both 
the lowest total stellar mass and the highest central angular mo- 
mentum. However, this trend does have significant scatter since 
other processes such as turbulent angular momentum transport 
and /V-body dynamics will affect the growth rate of the stellar 
cluster. 

4.3. Evolution of Protostellar Disk 

As the gas collapses and begins to form sink particles, it 
also develops into a flattened disk structure (Figure 9). We 
may estimate the dependence of the mass inflow rate inside the 
accretion disk on density and temperature using the following: 

M disk = 37rv2, (16) 

where £ is the disk surface density and v is estimated based upon 
the prescription introduced by Shakura & Sunyaev (1973), 

v = a S sH v c s . (17) 

// p is the pressure scale height of the disk and ass is a 
non-dimensional parameter ranging between ~10 -2 and 1, 
depending on the nature of angular momentum transport in the 
disk. In Figure 10 we show the inner surface density profile as 
measured within a thin 20 AU slice through the central disk, 
measured just prior to the formation of the first sink particle. 
We compare with the same set of other minihalos as previously 
discussed. Note that the length of 20 AU is chosen because this is 
the resolution limit of the comparison minihalos first presented 
in Stacy & Bromm (2013). Even on scales as small as 100 AU, 
the disk has surface densities as much as three times lower than 
the comparison minihalos. From Equation (16) we see that this 
leads to a correspondingly reduced disk accretion rate which 
persists well after the first sink appears. 

Let us roughly estimate ass — 0.1 and c s ~ 3 km s -1 , 
corresponding to a temperature of ~1000 K. The disk scale 
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Figure 7. Left: enclosed mass M en c vs. radius for various minihalos at the point of first sink or protostar formation. Right: ratio of M en c to Bonnor-Ebert mass Mbe 
vs. M e nc - The solid red line represents the simulation presented in this work just prior to initial sink formation, while the dashed red line represents this simulation 
5000 yr after sink formation. Other lines have the same meanings as in previous figures. The central point is taken as the most dense gas particle or most massive sink. 
Minihalos with higher overall accretion rates and less rotational support have lower overall M enc at a given radius and a lower ratio of M enc to A/be- 
(A color version of this figure is available in the online journal.) 
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Figure 8. Total mass accreted onto the stellar system after 5000 yr vs. the 
total angular momentum within the central 200 M© of enclosed mass. Angular 
momentum is measured just before the first sink forms. Values are taken from 
this simulation (red asterisk) as well as the minihalos studied in Stacy & Bromm 
(2013; filled circles). Note that the gas with the highest angular momentum also 
has the lowest stellar mass and accretion rate. 

(A color version of this figure is available in the online journal.) 



height at R — 100 AU may be estimated as H p /R ~ c s /v nA { R). 
From Figure 2 we estimate v ml (R =100 AU) ~ 3 km s -1 and 
thus H p ~ 100 AU. We then find v ~ 5 x 10 19 cm 2 s _1 . 

The surface density at 100 AU ranges from 30 to 100 g cm -1 , 
yielding Mdisk rates which range from 2 x 1 0 4 to 6 x 
10~ 4 M q yr -1 . Assuming the sinks grow at similar rates, we 
would expect after 10,000 yr that the stellar systems would 
reach a total mass of 2-6 M G . The lower end of this mass range 
is within good agreement with the total mass accretion rate seen 
in the simulation. However, the upper end is still somewhat 
lower than the total sink masses found in other calculations 
(Figure 8). Further variations in Mdisk are likely to come from 
differences in c s (i.e., warmer gas temperatures) and ass- This 
also indicates the approximative nature of our calculation. 


The range of accretion rates between the stellar disks also 
widens over time, such that after 5000 yr the total disk mass 
of the fastest-accreting halo from Stacy & Bromm (2013) is 
a factor of ten greater than the total disk mass within our z = 
1 5 minihalo. To illustrate this, in Figure 1 1 we show the resulting 
evolution of the disk mass, as well as the mass of disks taken 
from the simulations of Stacy & Bromm (2013) and Greif et al. 
(2012). In this figure, time is measured with respect to when 
the first sink forms. Note that the disk evolution of the green 
lines does not extend beyond a time of zero because these were 
taken from simulations which did not form sinks. The particular 
point that marks the transition from the stellar disk to the outer 
envelope is somewhat ambiguous, so to determine whether a 
gas particle is part of the disk we choose the simple criterion 
that it must have n > 10 9 cm 3 and /h 2 > 10 3 . Therefore, 
only dense and molecular gas is included. Compared to other 
studies, the z = 1 5 stellar disk grows at a very low rate. The slow 
overall growth rate of the mass of minihalo gas has translated to a 
slowly growing disk as well. The second-most slowly growing 
disk (green-dotted line in Figure 11) similarly belongs to the 
second-most slowly growing comparison minihalo. 

We briefly note that as mass falls into the sink, a portion of the 
gas gets heated to the virial temperature 7’ vlr of the sink through 
release of gravitational potential (Figure 1). Given a sink mass 
of 0.2 M g and using the accretion radius of 1.0 AU, we find a 
temperature of 


T vir ~ GMsinkmH ~ io 4 K. 

acc 


This “warm bubble” of neutral gas expands at its sound speed 
of c s < 10 km s _1 . Thus, by ~5000 yr the warm bubble has 
reached a distances of approximately c s t = 10,000 AU, which 
corresponds to gas of density n ~ 10 7 cm -3 . This warm phase 
of gas is visible in Figure 1 (see Turk et al. 2010 for further 
discussion of this warm and neutral gas phase). 

These disk properties lead to very low sink accretion rates 
(Figure 12). The disk gas fragments to form a second and third 
sink 600 and 900 yr after the first sink has formed, and two 
more between 3000 and 4000 yr. During the simulation two 


10 


The Astrophysical Journal, 785:73 (18pp), 2014 April 10 


Stacy & Bromm 



Size: 500 AU 



Size: 500 AU 


Figure 9. Density projection of central 500 AU of gas around first sink. Top and bottom rows are the x—z and x—y planes, respectively. From left to right, times after 
sink formation are 3, 2000, and 3500 yr. The asterisk marks the largest sink. The plus symbol marks the second largest sink. The diamonds are other secondary sinks. 
Note the rapid changes in the sink orbital motion and the structure of the protostellar disk. 

(A color version of this figure is available in the online journal.) 



Figure 10. Surface density profile of this simulation (solid red line) as well as 
two selected minihalos from Stacy & Bromm (2013; black dotted and dashed 
lines) and two from Greif et al. (2012; green lines). Profiles are shown just prior 
to the initial sink formation. Note that the Stacy & Bromm (2013) calculations 
were resolved down to ~10 AU, so we cannot calculate their surface densities 
on scales smaller than this. 

(A color version of this figure is available in the online journal.) 

other sinks form but quickly merge with the initial sink. A least- 
squares power-low fit to the growth rate of the three largest sinks 
remaining at f acc ~ 5000 yr yields 

M*,! ~ 0.61 M o (f/1000yr)°' 23 (18) 


0.51 M Q (t / 1000 yr) 028 

(19) 

0.62 M Q (t / 1000 yr) 032 

(20) 


Extending these power laws to 1 Myr, a typical accretion time 
for low-mass stars, these sinks would then reach 3.1, 3.5, and 
5.7 M q . A similar approximation for the growth of the least- 
massive sink predicts a late-time mass of 0.4 M 0 . 

It is uncertain how many more protostars will form at later 
times and what masses they would reach, since we do not 
follow the simulation for sufficient time to track the longer- 
term evolution of the disk and surrounding ~1000M e core as 
the protostars grow. As can be seen in Figure 7, the ratio of 
M enc to Mbe has a peak at M enc < 1000 Mq and then a rapid 
drop-off on larger scales. A similar drop off is seen at greater 
M enc for the more rapidly accreting halos. This central several 
hundred solar masses of material gravitationally infalls toward 
the central regions at a relatively slow rate of < 1 (T 3 M 0 yr -1 , 
leading to little change in the outer M enc profile after 5000 yr 
(see red dashed lines in Figure 7). 

At the same time, the protostars followed in our simulation 
are not projected to become massive enough to develop an H ii 
region that will blow away the gas. However, the luminosity 
and LW emission of the protostars will still serve to heat and 
stabilize the central gas, and the warm phase of the dense gas will 
continue to grow (Figure 1). We define the dense warm phase 
as gas which has n > 10 9 cnT 3 and /h 2 < 10~ 3 , such that 
only non-molecular gas is included. This phase approximately 
consists of ~2 M 0 of gas by the end of the simulation. The 
diversion of gas to the warm phase instead of the cool disk 
also helps to explain why M disk does not continue to increase 
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Figure 11. Evolution of disk mass over time. Time is measured with respect to 
the point at which a sink or protostar first forms. The solid red line represents the 
disk mass for our z = 15 minihalo. The black dotted line represents the minihalo 
from Stacy & Bromm (2013) that had the lowest overall stellar accretion rate, 
while the black dashed line represents that which had the greatest accretion 
rate. Green dotted and dashed lines additionally show the disk mass within two 
minihalos from Greif et al. (2012) which had the lowest and highest accretion 
rates, respectively. 

(A color version of this figure is available in the online journal.) 



above ~20 Mq after ~2000 yr of sink accretion. Even with such 
feedback effects, however, we cannot rule out the possibility 
that over the subsequent 10 5 to 10 6 yr, gas inflow onto the 
disk will continue until a star reaches ~10 Mq and ionizes its 
surroundings. We conjecture that stellar masses greater than 
10 Mq will not be necessary to halt slow inflow like that seen 
in our simulation, and that even if we followed our simulation 
for very long times we still would not see a Pop III star reach 
the more typical masses of 50-100 Mq. This will be confirmed 
with future numerical work. 

4.4. Evolution of Stellar Orbits 

During the disk evolution, the distance of the secondary sinks 
from the main sink ranges from ~ 1 AU to a few hundred AU. 


In comparison, the occurrence of Roche-lobe overflow requires 
the size of one of the binary members to exceed its Roche-lobe 
radius r t ;. 


— = max 
a 


0.46224 


\ 1/3 

. ) , 0.38 + 0.2 log 10 g 

i +qJ 


(21) 


for 0 < q < 0.8 (Paczynski 1971), where q is the binary 
mass ratio, and a is the semi-major axis. For q — 0. 1 we have 
ri = 0.2a, while q = 1 yields >q = 0.38a. This brackets 
the expected range of mass ratios for our simulated binary 
system. For a as small as 1 AU, we have ri ~ 3-6 x 10 12 cm, 
or 40-80 Rq. As AGB stars can reach well over 1 00 R (:) , this 
highlights the possibility that an AGB Pop III star may transfer 
mass to its companion. 

This is a particularly interesting possibility for the lowest- 
mass star, since it may experience additional close encounters 
as it orbits through and around the stellar system. If such a star 
remains below the “survival threshold” of 0.8 Mq (Figure 12), 
it may be observable in the present-day as a primordial AGB- 
enriched star in the Milky Way halo or nearby dwarf galaxy. This 
smallest sink experiences close encounters with the largest sink 
at,e.g., 1500 and 2500 yr (red line in the right panel of Figure 12) 
that nearly eject it from the disk. These encounters occur when 
the sink has grown to only 0.25 M Q and has a velocity of 
~5 km s" 1 relative to the disk. This is not quite sufficient to 
escape the stellar system, however. It is uncertain how much 
more it will grow as it continues its orbit through the accretion 
disk, but it is still only 0.25 M 0 at the end of the simulation and 
may remain below 1 M Q over its main-sequence lifetime (see 
also e.g., Johnson & Khochfar 201 1 for further discussion). 

We emphasize the speculative nature of the above scenario. 
An AGB phase for the larger protostars of our system would 
not occur until ~10 8 yr. It is at this later time that the smallest 
protostar’s orbit would need to come within > 1 AU for mass 
transfer to occur. While its orbit over the first 5000 yr ranges 
between one and a few hundred AU, it remains uncertain how the 
orbit will evolve over the much longer AGB timescales, whether 
the star will undergo an ejection before this time, etc. However, 
in our simulation we do still see the basic initial requirements 
for our scenario: that a slowly accreting low-mass star is in close 



time [yr] time [yr] 

Figure 12. Left: sink growth over time. The solid line represents the first and largest sink. The dotted line is the growth of the second-largest sink, while the other 
three black lines depict the growth of the three remaining sinks that survive to the end of our simulation. The red line depicts the total sink mass over time. The blue 
dash-dotted line depicts the “survival threshold,” the maximum mass for a star that could survive to the present-day. Right: distance of secondary sinks from the most 
massive sink over time. Line styles refer to the same sinks as in the left panel, but with different colors for more visible contrast between lines. 

(A color version of this figure is available in the online journal.) 
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Figure 13. Evolution of various properties of near-sink gas: Left: accretion rate onto the main sink. Right: value of a, where a = 1 corresponds to gas flowing radially 
toward the sink, while a = 0 corresponds to a rotationally dominated flow. The sink initially forms from an already disk-like gas configuration, so a is initially very 
low. Prior to 2500 yr, phases of more spherical motion around the sink generally correspond to periods of more rapid accretion. 


orbit around a larger star on track to eventually undergo an AGB 
phase. 

The variable orbital motion of the sinks further contributes 
to the high variability of the sink accretion rates. The accretion 
rate onto the main sink, shown in Figure 13, is nearly M ~ 
10 -2 M e yr -1 for the first few hundred years, but quickly 
drops to <10~ 3 Mq yr -1 with periods where M ~ 0. As the 
sink orbits through the stellar disk, the value of a is similarly 
variable, where a = 1 corresponds to radially dominated gas 
motion toward the sink and a = 0 corresponds to a rotationally 
dominated flow. Prior to 2500 yr, periods where a is closer to 
one corresponds to periods of more rapid accretion. In the latter 
half of the simulation when further disk fragmentation occurs, 
a becomes significantly more variable. On average, the main 
sink accretes at ~2 x 10~ 4 M G yr -1 ; 

5. INFLUENCE OF THE LYMAN-WERNER 
BACKGROUND 

5.1. Overview 

A photodisssociating LW background built up by earlier- 
forming Pop III stars may slow or prevent the cooling and 
collapse of the gas in our z — 15 minihalo. However, the effect 
of such a background remains very uncertain. We first note that 
our analysis shows it is not simply the redshift but also the 
rotational structure of the gas that drives the unusually slow 
infall rate onto this minihalo. There is substantial variation in 
minihalo characteristics seen at all redshifts. This implies that 
such a highly rotationally supported, slowly accreting gas cloud 
may also exist within some z = 20-30 minihalos. 

When considering the effect of LW radiation, it is indeed 
appropriate to focus on the global background radiation, as 
opposed to radiation from a particular source. Given our box 
size, the nearest minihalo “outside” of the box would be 200 kpc 
(comoving) away, or ~20 kpc (physical) away. Works by, e.g., 
Dijkstra et al. (2008) and Johnson et al. (2013) find that it 
is only rare high-density regions where LW flux from local 
sources will dominate over the global background, not regions 
like that in our simulation. However, the quickness with which 
this background will grow is very uncertain, and depends upon 
the early Pop III IMF For instance, the semi-analytic models 
of Crosby et al. (2013) find that J21 is at least ~1 at z = 15, 
where /21 represents units of 10~ 21 erg s” 1 cnT 2 Hz -1 sr _1 . The 
simulations of Johnson et al. (2013), on the other hand, find at 



Figure 14. Evolution of /shield with redshift within the gas core of the minihalo. 
Dotted line also shows the maximum gas density « max . 


this redshift that the overall J 21 is an order of magnitude lower, 
/21 ~0.1. 

It is furthermore possible that the concurrently growing 
X-ray background from Pop III remnant black holes could have 
an opposing effect to the LW background. For instance, Jeon 
et al. (2012) find that the X-ray radiation from a high-mass BH 
binary can provide positive feedback such that gas collapse into 
distant minihalos is facilitated via Hi cooling promoted by the 
strong X-ray emission. In contrast, semianalytic models by, e.g., 
Tanaka et al. (2012) find that the X-ray heating of the IGM is a 
stronger effect than the associated enhancement of H 2 formation 
in potential star-forming regions. 

Simulations such as those in Machacek et al. (2003) and 
O’Shea & Norman (2008), which found that the LW background 
delays gas collapse, did not account for the effect of H 2 self- 
shielding. They furthermore assumed a /21 that was constant 
instead of gradually evolving. To properly include the effects of 
this background, we would need to understand the still uncertain 
rate at which these backgrounds build up, and we also would 
need to account for H 2 self-shielding within the minihalo. 

The value of /shield represents the factor by which H 2 ab- 
sorption from the IGM and the outer parts of the minihalo will 
reduce the local LW flux within the inner star-forming parts of 
the minihalo. Figure 14 shows estimates of the average /shield 
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within the gas core of the minihalo, defined as the central gas 
with densities within a factor of ten of the maximum gas den- 
sity H m ax- We estimate /shield in the same manner as in Johnson 
et al. (2013), in turn based upon Draine & Bertoldi (1996) and 
Wolcott-Green et al. (2011): 


/shield (Wt 2 , T ) 


0.965 0.035 

(1 + x/b^) 11 (l+.r) 0 - 5 

x exp[— 8.5 x 10 4 (1 +x) 0 ' 5 ], (22) 


where x = 1Vh 2 /5x 10 14 cm 4 , b$ = b/ 10 5 cm s , and b is the 
Doppler broadening parameter given by b = ik^T/m^) 1 ^ 1 (see 
Johnson et al. 2013 for further details). By z < 25, /shield will 
reduce the local LW flux by over an order of magnitude, greatly 
helping to reduce any possible effects of the LW background. 

We finally note that we do not argue that Pop III stars were 
typically low mass, but that in rare environments, such a low- 
mass formation mode could occur. This is potentially important 
as it would allow observers to detect such Pop III fossils 
as surviving stars in our local cosmic neighborhood (stellar 
archaeology). The standard picture to date has been that Pop III 
stars always grow to masses that would have led to their death 
a long time ago. Because the strength of the LW background is 
subject to huge uncertainties, including spatial fluctuations and 
local opacity, conditions such as those simulated here cannot be 
excluded. 


5.2. Numerical Tests 

We further numerically examine LW effects with a set 
of re-simulations of the initial minihalo collapse, beginning 
from z = 50, but this time including a LW background. 
We follow the gas collapse up to the point just before sink 
formation ( n = 10 16 cm' 3 ), but we did not have sufficient 
computational resources to follow the gas evolution further. The 
LW background grows with time as 

J 21 = Ju,o x Kr ( ^» )/5 , (23) 

where J 2 1,0 is a normalization parameter set to range from 0.1 
to 10, and zo = 10. When J 210 = 1, we obtain a good fit to 
the LW background evolution presented in Figure 1 of Greif & 
Bromm (2006; see also Pawlik et al. 2013). From this we apply 
a photo-dissociation rate of 

k Hl = 1.38 x 10- 12 /shield J 21 s _1 . (24) 

We find that the LW background initially serves to delay 
minihalo collapse. For the more extreme 721,0 = 10 case, 
the minihalo gas still has not collapsed at z = 13.3. At this 
time the densest gas has densities of only n = 10 cm -3 . To 
save computational time we did not follow this simulation 
further. However, we conclude that at such high 721,0 values the 
gas evolution is significantly altered. Gas collapse will likely 
occur only once the minihalo reaches higher mass and virial 
temperatures of ~10 4 K. Gas fragmentation may be suppressed, 
and the gas may collapse directly into a black hole (see, e.g.. Oh 
& Haiman 2002; Bromm & Loeb 2003). 

The gas evolution is much less affected for the more physi- 
cally realistic value of 721,0 = 0. 1, in which case collapse to sink 
particle densities of n = 10 16 cm -3 is delayed by 1.7 x 10 6 yr. 
In Figure 15 we show the state of the gas at the point that the 
gas has reached n = 10 16 cm -3 in each test simulation. For 
721,0 = 0.1 (solid black line in Figure 15), the gas profile does 


not significantly differ from the fiducial 721,0 = 0 case (red line 
in Figure 15). As expected, outside of 10 5 AU the H 2 fraction 
is reduced due to the LW background, while shielding is effec- 
tive in the more central regions and even allows for a slightly 
larger H 2 fraction inside of 10 5 AU. In the inner 1000 AU, 
the gas properties do not differ by more than a few tens of 
percent, though the central densities and estimated spherical ac- 
cretion rate (see Equation (14)) are somewhat enhanced. The 
enclosed mass at all given radii, as well as the ratio of M enc to 
M b e, are also slightly larger (Figure 16). This leads to a slightly 
higher overall disk mass as the gas approaches n = 10 16 cm -3 , 
when Mdisk ~ 20 M 0 as opposed to ~ 1 6 M 0 in the fiducial case 
(Figure 1 7). Note, however, that these values of M^k and M sp he re 
are approximately half to a tenth as large as those found in the 
other minihalos discussed in Section 4.2. Thus, under a small 
721,0 = 0.1 background we would still likely find an unusually 
low-mass Pop III system, though probably more massive than 
the fiducial case. 

For 72i,o = 1, the differences are more significant (dashed 
line in Figure 15). Collapse is delayed by 3.3 x 10 7 yr before 
the gas finally reaches the sink density of n — 10 16 cm '. The 
H 2 fraction is approximately half of that found in the fiducial 
case at distances greater than 10 4 AU. However, inside 10 4 AU 
shielding becomes very effective, and the H 2 increases more 
rapidly than the other test cases as the radius declines. In the 
central 10 4 AU the gas density, sound speed, and infall velocity 
are all reduced, and Msphere is up to an order of magnitude 
below the fiducial case. The enclosed mass at all given radii 
is also reduced by up to a factor of a few (Figure 16). The 
central 30 M enc has a higher ratio of M enc to M B e, but outside 
of this region the ratio is much reduced and the gas is more 
gravitationally stable. Overall, this leads to M^sk ~ 5 M 0 when 
n = 10 16 cm 3 , several times smaller than the other cases. 

When comparing these simulations at the same peak densities, 
but not at the same physical time, the effect of increasing the 
LW background is not monotonic. For the more realistic values 
of ./ 2 i,o = 0. 1 or 1, we find central accretion rates that are either 
slightly enhanced or significantly reduced. In these cases the 
LW background does not seem to change our general finding 
of a minihalo which hosts a system of very slowly accreting 
Pop III stars. 

6. DISCUSSION AND CONCLUSIONS 

We present a three-dimensional simulation of the formation 
and growth of a Pop III stellar system. This calculation was 
initialized on cosmological scales while resolving lengths as 
small as 1 AU. We found that the host minihalo formed at 
an unusually low redshift of z — 15, leading to a low DM and 
baryonic accretion rate as well as a low-mass and slowly growing 
stellar system. The stars in the system can be expected to reach 
~0.5 to 5 Mq after 1 Myr. This is nearly an order of magnitude 
slower than rates typically found within z > 20 stellar systems 
(e.g., Stacy et al. 2010; Greif et al. 2011; Smith et al. 2011). 
We additionally find that a LW background as high as 721,0 = 1 
will not significantly change our finding of uncharacteristically 
small infall rates onto the central star-forming region of our 
minihalo. 

It is uncertain how common such low-mass Pop III systems 
will be. The minihalo we present here hosts the slowest- 
accreting Pop III system of approximately 10 minihalos which 
are included in our comparison in Section 4. Hirano et al. 
(2014) use approximately 100 minihalos from cosmological 
simulations to initialize two-dimensional simulations of Pop III 
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Figure 15. Radial profiles of the central gas within our simulated minihalo under various strengths of LW background, measured when the gas has reached a density 
of n = 10 16 cm -3 in the respective simulations. The solid black lines denote the 721,0 = 0.1 case, while dashed lines represent 721,0 = 1. The red lines are taken from 
the fiducial 72i,o = 0 simulation discussed throughout this work. Top left: H 2 fraction. Top right: electron fraction. Middle left: number density. Middle right: sound 
speed. Bottom left: radial infall velocity. Bottom right: estimated spherical accretion rate (M S phe re ; see Equation (14)). The gas profiles are not significantly altered for 
721,0 = 0.1, though density, infall velocity, and M sp here are slightly enhanced. The 721,0 = 1 background surprisingly leads to the generally opposite effect of reduced 
densities, sound speeds, and infall rates. 

(A color version of this figure is available in the online journal.) 


stellar growth under feedback. Only one star can be followed 
per minihalo, and in this case their smallest star is expected to 
grow to ~9 Mq. Were they able to follow fragmentation within 
this particular minihalo, it is possible that the stellar mass may 
have instead been distributed among several stars of lower mass. 
We thus make a very rough estimate that one out of a few tens 
to one out of 100 minihalos will host a low-mass system similar 
to what we find in our simulation. Considering that the mass 


of 10 5 — 10 6 minihalos will ultimately become incorporated into 
a Milky Way-type galaxy, it is conceivable that on the order 
of thousands of Pop III stars from such low-mass systems may 
exist in the nearby Galactic halo. 

Further study will be necessary to determine more precisely 
how common such low-mass Pop III systems are and whether 
our 2=15 Pop III system indicates a more general transition 
in the Pop III IMF from z ~ 30 to ~ 10, or if our particular 
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Figure 16. Top: enclosed mass M enc vs. radius for our test cases at the point 
when the maximum gas density first reaches 10 16 cm -3 . The solid black lines 
denote the J 2 i,o = 0.1 case, while dashed lines represent 721,0 = 1. The red 
lines are taken from the fiducial 721,0 = 0 simulation discussed throughout this 
work. The central point is taken to be the most dense gas particle. Bottom: ratio 
of Afenc to A7 be vs. M en c- Lines have same meaning as in the upper panel. While 
the 72i,o = 0.1 case has slightly enhanced mass and M enc to Mbe ratio, these 
are somewhat reduced for the 72i,o = 1 case. 

(A color version of this figure is available in the online journal.) 


minihalo was an unusual case even for its low redshift. While 
minihalo comparisons of other works do not find this same 
transition in accretion rate (Gao et al. 2007; O’Shea & Norman 
2007), the numerical analyses by de Souza et al. (2013) do 
find a transition in the spin distribution of gas within minihalos 
as redshift declines. According to the semi-analytic model of 
McKee & Tan (2008), such differences in rotational support 
will in turn lead to differences in the protostellar feedback and 
the final mass reached by the protostar, de Souza et al. (2013) 
argue that this will cause the Pop III IMF peak to shift to lower 
mass with lower redshift. This remains to be tested with more 
physically detailed simulations. 

Understanding the Pop III IMF evolution with redshift will 
additionally require a greater knowledge of how the build- 
up of global background radiation as well as magnetic fields 
proceeded at this redshift (e.g., Schleicher et al. 2010; Schober 
et al. 2012; Turk et al. 2012). Understanding how the IMF 
evolved to later times is particularly important given that recent 



Figure 17. Evolution of disk mass over time for 721,0 = 0, 721,0 = 0.1, and 
721,0 = L Line styles have the same meaning as in the previous figure. Time is 
measured with respect to the time when the gas first reaches n = 10 16 cm -3 . 
The 72i,o = 0.1 background leads to a slightly enhanced M^k as compared 
with the fiducial 721,0 = 0 model. The larger 721,0 = 1 background in fact 
leads to Mdisk which is consistently several times smaller than the fiducial case. 
Realistic LW backgrounds are therefore unlikely to change our overall finding 
of unusually low Pop III accretion rates within our simulated minihalo. 

(A color version of this figure is available in the online journal.) 


work has indicated that metal-free gas will indeed survive 
to relatively low redshift. For instance, Simcoe et al. (2012) 
reported observations of extremely low metallicity or possibly 
metal-free gas within a z ~ 7 damped Lya system, while 
Fumagalli et al. (2011) reported the detection of metal-free 
gas within Lyman-limit systems at z > 3. Numerical work 
(e.g., Muratov et al. 2013; see also Scannapieco et al. 2005) has 
similarly indicated that Pop III star formation can continue to 
Z ~ 6. Our preliminary tests presented in Section 5, however, 
indicate that for a range of LW backgrounds our Pop III system 
will still undergo unusually low accretion rates. 

AGB stars are important to understanding the early chemi- 
cal evolution of the galaxy (e.g., Karakas 2010; Campbell & 
Lattanzio 2008; Karakas & Lugaro 2010). They are known 
to produce significant quantities of carbon, nitrogen, and 
s-process elements. (Busso et al. 2001; Siess et al. 2002; Siess & 
Goriely 2003). They may also help to explain the observation of 
large amounts of dust in high-/ galaxies and quasars, which im- 
plies rapid dust production in the early universe (Bertoldi et al. 
2003; Valiante et al. 2009; Gall et al. 201 1). Our results reveal a 
pathway for such AGB star formation within primordial gas. In 
particular, the larger 3-5 M 0 stars of our simulation are suffi- 
ciently short-lived to undergo a metal and dust-enriching AGB 
phase by z > 6, while the smaller < 1 M Q star is long-lived 
enough that it may still be observed as carrying the enrichment 
signatures of its larger companions. The physical scenario sug- 
gested by our simulations may indeed explain the abundances 
of certain metal-poor stars in the Milky Way halo, particularly 
those with unusual features in C, N, and O. Studies have found 
that some of these may be Pop III stars which received material 
from an AGB companion that is now a white dwarf (e.g., Suda 
et al. 2004, 2013). 

Studies over recent years have found increasing complexity 
in the nature of Pop III stars. Though dust and metallicity 
by definition do not play a role in primordial star formation, 
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other physical processes such as multiplicity (e.g., Clark et al. 
2008; Stacy et al. 2010; Greif et al. 2011) as well as the 
binary nature and rotation rate of Pop III stars (e.g., Stacy & 
Bromm 2013) will also be of crucial importance. In addition, 
feedback (e.g., Hosokawa et al. 2011; Smith et al. 2011; Stacy 
et al. 2012) and magnetic fields play a central role in Pop III 
protostellar accretion. However, our work shows that even when 
such processes are not included in simulations, the differences 
in minihalo environments alone can lead to substantial variation 
between primordial stellar clusters. As shown by Jappsen et al. 
(2009; see also Dopcke et al. 2013), our work demonstrates that 
the transition to a low-mass IMF will depend upon not only a 
critical metallicity (e.g., Bromm et al. 2001), but also upon the 
characteristics of the primordial star-forming region. 

We thus predict the rare existence of low-mass Pop III stars 
that have survived until the present day and that may show 
evidence of enrichment from a companion’s AGB-phase mass 
overflow. Continued improvements in simulation techniques and 
computational power, as well as constraints provided by obser- 
vations such as abundance measurements of metal-poor stars in 
the Galactic halo and dwarf galaxies (e.g., Beers & Christlieb 
2005; Frebel et al. 2005; Caffau et al. 2011; see also Karlsson 
et al. 2013), will provide an increasingly refined picture of the 
role Pop III stars played in shaping the early universe. Ulti- 
mately, the question of whether true Pop III survivors exist, e.g., 
in the guise of AGB self-enrichment as suggested in this paper, is 
a question that can be tested empirically. The potential for such 
stellar archaeological constraints is demonstrated by the excit- 
ing recent discovery by the SkyMapper Southern Sky Survey of 
a star with no detected Fe-peak elements, but low abundances of 
C, N, and O (Keller et al. 2014). Indeed, the detected abundance 
pattern closely resembles an AGB self-enrichment pattern with 
the one crucial exception of an extremely low, but non-zero, 
Ca abundance. The latter cannot be accommodated with AGB 
enrichment, but instead points to the signature of supernova 
enrichment from a massive Pop III progenitor, rendering this 
extreme star second-generation. However, the search for Pop III 
fossils is clearly within the reach of current and upcoming sur- 
veys. 

The authors thank the anonymous referee who helped them 
improve this manuscript. A.S. is grateful for support from the 
JWST Postdoctoral Fellowship through the NASA Postdoc- 
toral Program (NPP). V.B. acknowledges support from NASA 
through Astrophysics Theory and Fundamental Physics Pro- 
gram grant NNX09AJ33G and from NSF through grant AST- 
1009928. Resources supporting this work were provided by 
the NASA High-End Computing (HEC) Program through the 
NASA Advanced Supercomputing (NAS) Division at Ames 
Research Center. 

REFERENCES 

Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, NewA, 2, 181 
Abel, T„ Bryan, G. L., & Norman, M. L. 2002, Sci, 295, 93 
Alvarez, M. A., Bromm, V., & Shapiro, P. R. 2006, ApJ, 639, 621 
Bate, M. R„ Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362 
Bate, M. R., & Burkert, A. 1997, MNRAS, 288, 1060 
Beers, T. C., & Christlieb, N. 2005, ARA&A, 43, 531 
Bertoldi, F„ Carilli, C. L.. Cox, P, et al. 2003, A&A, 406, L55 
Bromm, V. 2013, RPPh, 76, 112901 

Bromm, V., Coppi, P. S., & Larson, R. B. 2002, ApJ, 564, 23 

Bromm, V., Ferrara, A., Coppi, P. S., & Larson, R. B. 2001, MNRAS, 328, 969 

Bromm, V., & Loeb, A. 2003, ApJ, 596, 34 

Bromm, V., & Loeb, A. 2004, NewA, 9, 353 

Bromm, V., Yoshida, N., & Hemquist, L. 2003, ApJL, 596, L135 


Busso, M., Gallino, R., Lambert, D. L., Travaglio, C., & Smith, V. V. 2001, ApJ, 
557, 802 

Caffau, E., et al. 201 1, Natur, 477, 67 
Campbell, S. W„ & Lattanzio, J. C. 2008, A&A, 490, 769 
Chatzopoulos, E., & Wheeler, J. C. 2012, ApJ, 748, 42 
Cherchneff, I., & Dwek, E. 2010, ApJ, 713, 1 
Clark, P. C., Glover, S. C. O., & Klessen, R. S. 2008, ApJ, 672, 757 
Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bromm, V. 2011a, ApJ, 
727. 110 

Clark, P. C., Glover, S. C. O., Smith, R. J., et al. 2011b, Sci, 331, 1040 
Crosby, B. D„ O’Shea, B. W„ Smith, B. D„ Turk, M. J., & Hahn, O. 2013, ApJ, 
773, 108 

Dekel, A., Birnboim, Y„ Engel, G„ et al. 2009, Natur, 457, 45 1 
de Souza, R. S., Ciardi, B., Maio, U., & Ferrara, A. 2013, MNRAS, 428, 2109 
Dijkstra, M., Haiman, Z., Mesinger, A., & Wyithe, J. S. B. 2008, MNRAS, 
391. 1961 

Dopcke, G., Glover, S. C. O., Clark, P. C., & Klessen, R. S. 2013, ApJ, 766, 103 
Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269 

Faucher-Giguere, C.-A., Keres, D„ & Ma, C.-P. 2011, MNRAS, 417, 2982 
Frebel, A., Aoki, W., Christlieb, N.. et al. 2005, Natur, 434, 871 
Frommhold, L. 1994, Collision-induced Absorption in Gases (Cambridge: 
Cambridge Univ. Press) 

Fumagalli, M., O'Meara, J. M., & Prochaska, J. X. 2011, Sci, 334, 1245 
Gall, C., Hjorth, J., & Andersen, A. C. 2011, A&ARv, 19, 43 
Gao, L„ White, S. D. M„ Jenkins, A., Frenk, C. S., & Springel, V. 2005, MNRAS, 
363, 379 

Gao, L„ Yoshida, N„ Abel, T„ et al. 2007, MNRAS, 378, 449 
Genel, S., Genzel, R., Bouche, N., et al. 2008, ApJ, 688, 789 
Greif, T„ Springel, V., White, S., et al. 2011, ApJ, 737, 75 
Greif, T. H., & Bromm, V. 2006, MNRAS, 373, 128 
Greif, T. H„ Bromm, V., Clark, P. C., et al. 2012, MNRAS, 424, 399 
Greif, T. H„ Glover, S. C. O., Bromm, V., & Klessen, R. S. 2010, ApJ, 716, 510 
Greif, T. H., Johnson, J. L., Bromm, V., & Klessen, R. S. 2007, ApJ, 670, 1 
Greif, T. H., Johnson, J. L., Klessen, R. S., & Bromm, V. 2009, MNRAS, 
399, 639 

Greif, T. H., Springel, V., & Bromm, V. 2013, MNRAS, 434, 3408 
Haiman, Z., Thoul, A. A., & Loeb, A. 1996, ApJ, 464, 523 
Hansen, C. J., Kawaler, S. D., & Trimble, V. 2004, Stellar Interiors: Physical 
Principles, Structure, and Evolution (New York: Springer) 

Hartmann, L., Cassen, P, & Kenyon, S. J. 1997, ApJ, 475, 770 
Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532 
Henyey, L. G., Lelevier, R., & Levee, R. D. 1955, PASP, 67, 154 
Hirano, S., Hosokawa, T„ Yoshida, N., et al. 2014, ApJ, 781, 60 
Hirano, S., & Yoshida, N. 2013, ApJ, 763, 52 

Hosokawa, T„ Omukai, K„ Yoshida, N„ & Yorke, H. W. 2011, Sci, 334, 1250 
Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478 
Jappsen, A.-K., Klessen, R. S., Glover, S. C. O., & Mac Low, M.-M. 2009, ApJ, 
696, 1065 

Jeon, M„ Pawlik, A. H„ Greif, T. H„ et al. 2012, ApJ, 754, 34 

Johnson, J. L„ Dalla, V. C., & Khochfar, S. 2013, MNRAS, 428, 1857 

Johnson, J. L., Greif, T. H., & Bromm, V. 2007, ApJ, 665, 85 

Johnson, J. L„ & Khochfar, S. 2011, MNRAS, 413, 1184 

Karakas, A. I. 2010, MNRAS, 403, 1413 

Karakas, A. I., & Lugaro, M. 2010, PASA, 27, 227 

Karlsson, T., Bromm, V., & Bland-Hawthom, J. 2013, RvMP, 85, 809 

Kashlinsky, A., & Rees, M. J. 1983, MNRAS, 205, 955 

Keller, S. C„ Bessell, M. S„ Frebel, A., et al. 2014, Natur, 506, 463 

Kitayama, T„ Yoshida, N., Susa, H., & Umemura, M. 2004, ApJ, 613, 631 

Lacey, C., & Cole, S. 1993, MNRAS, 262, 627 

Latif, M. A., Schleicher, D. R. G., Schmidt, W., & Niemeyer, J. 2013, ApJL, 
772, L3 

Loeb, A. 2010, How Did the First Stars and Galaxies Form? (Princeton, NJ: 
Princeton Univ. Press) 

Machacek, M. E„ Bryan, G. L„ & Abel, T. 2003, MNRAS, 338, 273 

Madau, P., Ferrara, A., & Rees, M. J. 2001, ApJ, 555, 92 

Maio, U., Khochfar, S„ Johnson, J. L., & Ciardi, B. 2011, MNRAS, 414, 1 145 

Martel, H„ Evans, N. J„ & Shapiro, P. R. 2006, ApJS, 163, 122 

McKee, C. F„ & Tan, J. C. 2008, ApJ, 681, 771 

Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347 

Mori, M., Ferrara, A., & Madau, P. 2002, ApJ, 571, 40 

Muratov, A. L., Gnedin, O. Y., Gnedin, N. Y., & Zemp, M. 2013, ApJ, 773, 19 

Nakamura, F., & Umemura, M. 2001, ApJ, 548, 19 

Navarro, J. F„ & White, S. D. M. 1994, MNRAS, 267, 401 

Neistein, E., van den Bosch, F. C., & Dekel, A. 2006, MNRAS, 372, 933 

Nomoto, K„ Maeda, K„ Umeda, H., et al. 2003, PThPS, 151, 44 

Norman, M. L„ O'Shea, B. W., & Paschos, P. 2004, ApJL, 601, LI 15 

Oh, S. P., & Haiman, Z. 2002, ApJ, 569, 558 


17 


The Astrophysical Journal, 785:73 (18pp), 2014 April 10 


Stacy & Bromm 


Omukai, K„ & Palla, F. 2003, ApJ, 589, 677 
O’Shea, B. W„ & Norman, M. L. 2007, ApJ, 654, 66 
O’Shea, B. W„ & Norman, M. L. 2008, ApJ, 673, 14 
Paczynski, B. 1971, ARA&A, 9, 183 

Pawlik, A. H., Milosavljevic, M„ & Bromm, V. 2013, ApJ, 767, 59 

Press, W. H„ & Schechter, P. 1974, ApJ, 187, 425 

Prialnik, D„ & Livio, M. 1985, MNRAS, 216, 37 

Reed, D. S., Bower, R„ Frenk, C. S., et al. 2005, MNRAS, 363, 393 

Ripamonti, E„ & Abel, T. 2004, MNRAS, 348, 1019 

Ripamonti, E., Haardt, F., Ferrara, A., & Colpi, M. 2002, MNRAS, 
334, 401 

Scannapieco, C., Tissera, P. B., White, S. D. M., & Springel, V. 2005, MNRAS, 
364, 552 

Schleicher, D. R. G., Banerjee, R., Sur, S., et al. 2010, A&A, 522, Al 15 

Schober, J.. Schleicher, D., Federrath, C., et al. 2012, ApJ, 754, 99 

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 

Shu, F. H. 1977, ApJ, 214, 488 

Siess, L., & Goriely, S. 2003, NuPhA, 718, 524 

Siess, L., Livio, M., & Lattanzio, J. 2002, ApJ, 570, 329 

Simcoe, R. A., Sullivan. P. W., Cooksey, K. L., et al. 2012, Natur, 492, 79 

Smith, R. J., Glover, S. C. O., Clark, P. C., Greif, T., & Klessen, R. S. 

2011, MNRAS, 414, 3633 

Smith, R. J., Hosokawa, T., Omukai, K., Glover, S. C. O., & Klessen, R. S. 

2012, MNRAS, 424, 457 

Sokasian, A., Yoshida, N., Abel, T., Hernquist, L., & Springel, V. 2004, MNRAS, 
350, 47 

Springel, V. 2005, MNRAS, 364, 1 105 


Stacy, A., & Bromm, V. 2013, MNRAS, 433, 1094 
Stacy, A., Greif, T. H„ & Bromm, V. 2010, MNRAS, 403, 45 
Stacy, A., Greif, T. H„ & Bromm, V. 2012, MNRAS, 422, 290 
Stahler, S. W., Palla, F„ & Salpeter, E. E. 1986a, ApJ, 308, 697 
Stahler, S. W„ Palla, F., & Salpeter, E. E. 1986b, ApJ, 302, 590 
Suda, T., Aikawa, M., Machida, M. N., Fujimoto, M. Y., & Iben, I., Jr. 2004, ApJ, 
611,476 

Suda, T„ Komiya, Y„ Yamada, S., et al. 2013, MNRAS, 432, L46 
Tanaka, T„ Perna, R„ & Haiman, Z. 2012, MNRAS, 425, 2974 
Tegmark, M„ Silk, J., Rees, M. J., et al. 1997, ApJ, 474, 1 
Tortnen, G„ Bouchet, F. R„ & White, S. D. M. 1997, MNRAS, 286, 865 
Tomatore, L., Ferrara, A., & Schneider, R. 2007, MNRAS, 382, 945 
Turk, M. J., Abel, T„ & O'Shea, B. 2009, Sci, 325, 601 
Turk, M. J., Norman, M. L„ & Abel, T. 2010, ApJL, 725, L140 
Turk, M. J., Oishi, J. S., Abel, T„ & Bryan, G. L. 2012, ApJ, 745, 154 
Valiante, R., Schneider, R., Bianchi, S., & Andersen, A. C. 2009, MNRAS, 
397, 1661 

Wada, K„ & Venkatesan, A. 2003, ApJ, 591, 38 
Whalen, D„ Abel, T„ & Norman, M. L. 2004, ApJ, 610, 14 
Wise, J. H.. & Abel, T. 2008, ApJ, 685, 40 

Wise, J. H„ Turk, M. J., Norman, M. L„ & Abel, T. 2012, ApJ, 745, 50 
Wolcott-Green, J., & Haiman, Z. 201 1, MNRAS, 412, 2603 
Wolcott-Green, J., Haiman, Z., & Bryan, G. L. 201 1 , MNRAS, 418, 838 
Yoon, S.-C., Dierks, A., & Langer, N. 2012, A&A, 542, A113 
Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645 
Yoshida, N., Omukai, K., & Hernquist, L. 2008, Sci, 321, 669 
Yoshida, N., Omukai, K., Hernquist, L., & Abel, T. 2006, ApJ, 652, 6 


18 


